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The standard post-Newtonian approximation to gravitational waveforms, called T-approximants, 
from non-spinning black hole binaries are known not to be sufficiently accurate close to the last 
stable orbit of the system. A new approximation, called P-approximants, is believed to improve 
the accuracy of the waveforms rendering them applicable up to the last stable orbit. In this study 
we apply P-approximants to the case of a test-particle in equatorial orbit around a Kerr black hole 
parameterized by a spin parameter q that takes values between —1 and 1. In order to assess the 
performance of the two approximants we measure their effectualness (i.e. larger overlaps with the 
exact signal), and faithfulness (i.e. smaller biases while measuring the parameters of the signal) 
with the exact (numerical) waveforms. We find that in the case of prograde orbits, that is orbits 
whose angular momentum is in the same sense as the spin angular momentum of the black hole, 
T-approximant templates obtain an effectualness of ~ 0.99 for spins q < 0.75. For 0.75 < q < 0.95, 
the effectualness drops to about 0.82. The P-approximants achieve effectualness of > 0.99 for all 
spins up to q = 0.95. The bias in the estimation of parameters is much lower in the case of P- 
approximants than T-approximants. We find that P-approximants are both effectual and faithful 
and should be more effective than T-approximants as a detection template family when q > 0. For 
q < both T- and P-approximants perform equally well so that either of them could be used as a 
^ ' detection template family. 
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Stellar mass compact binaries consisting of double neutron stars (NS), double black holes (BH) or a mixed binary 
consisting of a neutron star and a black hole, are the primary targets for a direct first detection of gravitational waves 
(GW) by interferometric detectors, LIGO 0, VIRGO 0, GEO600 3], and TAMA 0|. Under radiation reaction the 

O^' orbit of a binary slowly decays, emitting a signal whose amplitude and frequency increases with time and is termed 
a "chirp" signal. While it is believed that there is a greater population of NS-NS binaries pa, l1 0, 0- 12 i it is the 

bi), BH-BH binaries that are the strongest candidates for detection since they can be seen from a greater volume, about 
two orders-of-magnitude greater than NS-NS binaries [a, [l(| . 
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I. INTRODUCTION 



A. Matched filtering 

In order to detect such sources one employs the method of matched filtering [ll| . Briefly, the method works as 
follows: Firstly, one creates a set of waveforms, or templates as they are called, that depend on a number of parameters 
of the source and its location and orientation relative to the detector. These templates are then cross-correlated with 
the detector output weighted by the inverse of the noise spectral density. If a signal, whose parameters are close to 
one of the template waveforms, is actually present in the detector output then the cross-correlation builds up, with 
the dominant contribution coming from frequencies where the noise spectral density is low. Thus, in the presence of 
a sufficiently strong signal the correlation will be much larger than the RMS correlation in the absence of any signal. 
How large should it be before we can be confident about the presence of a signal depends on the combination of the 
rate of inspiral events and the false alarm probability (see e.g. Ref. [Tj] for a simple estimation). 

The effectiveness of matched filtering depends on how well the phase evolution of the waveform is known. Even 
tiny instantaneous differences, as low as one part in 10 3 in the phase of the true signal that might be present in the 
detector output and the template that is used to dig it out could lead to a cumulative difference of several radians 
since one integrates over several hundreds to several thousands of cycles. In view of improving the signal-to-noise 
ratio for inspiral events there has been a world-wide effort in accurately computing the dynamics of a compact binary 
and the waveform it emits or to use phenomenologically defined detection template families [ljj, Hj, Il3 • 



B. Phasing of the coalescing binary signal 

There have been parallel efforts on using two different approximation schemes: On the one hand the post-Newtonian 
(PN) expansion of Einstein's equations has been used to treat the dynamics of two bodies of comparable masses with 
and without spin, in orbit around each other. This approximation is applicable when the velocities involved in the 
system are small but there is no restriction on the ratio of the masses [lg, UJl HH USi 1211 UlU UM ■ 

On the other hand, black hole perturbation theory has been used to compute the dynamics of a test particle in 
orbit around a spin-less or spinning black hole. Black hole perturbation theory does not make any assumptions on 
the velocity of the c omp one nts, but is valid only in the limit when the mass of one of the bodies is much less than 
the other jBLIM I2I 129. Wl\.H^ . 

The post-Newtonian approximation is a perturbative method which expands the equations of motion, binding 
energy and GW flux as a power series in v/c, where v is a typical velocity in the system and c is the speed of light. 
In the early stages of an inspiral, the radiation reaction time-scale trr ~ uj/Cj, where u> is the angular velocity and 
lu its time-derivative, is much greater than the orbital time-scale r or b ~ 1/lu. It is during this adiabatic regime that 
the post-Newtonian approximation works best. At present, the PN expansion for the case of comparable-masses is 
known to order O (w 6 ) [21| and O (u 7 ) |22| , for the energy and flux functions, respectively. However, at this order 
an arbitrary parameter exists in the expression for the flux. In order to see how well PN theory performs, we can 
compare two different systems. If we assume a NS-NS binary of masses (1.4, 1.4) M^ a nd a lower frequency cutoff 
of the detector at 40 Hz, then the "orbital velocity" of the binary is small, v ~ 0.12, [5_9j when it enters the detector 
bandwidth and the two stars are still largely separated, r ~ 70 M. The ratio of time-scales in the most sensitive 
regime of the detector is in the range 4.5 x 10 3 < TaR./T or b < 680. If on the other hand we take a BH-BH binary of 
masses (10, 10) M©, the orbital velocity is quite large, v ~ 0.23, and the separation is quite small, r ~ 19 M, upon 
entering the detector bandwidth. This is very close to the regime, v ~ 0.3, r ~ 11M, where the background curvature 
becomes strong and the motion relativistic. Once again, comparing time-scales, we obtain 140 < TRF>,/T or b < 40, where 
the final value is taken at the last stable orbit at /lso ~ 210 Hz. It is known that PN theory becomes inaccurate at 
an orbital separation of r < 10 M (29J. Therefore, post-Newtonian approximation becomes less valid for higher mass 
systems in the LIGO band but well describes the early stages of the inspiral of a NS-NS system visible in LIGO. 

As previously stated, black hole perturbation theory makes no assumptions about the orbital velocity of the compo- 
nents, but does restrict their masses. One assumes that a test particle of mass \x is in orbit about a central BH of mass 
M such that /i <C M. Assuming this restriction is satisfied we have an analytical expression for the energy. However, 
no analytical expression has been worked out for the gravitational wave flux emitted by such a system. Using black 
hole perturbation theory, a series approximation was initially calculated to O (v s ) by Poisson for a test particle in 
circular orbit around a Schwarzschild black hole [23]. The series was extended numerically to O (y 5 ) by Cutler et 
al. J2J], and then to O (u 8 ) by Tagoshi and Nakamura [2J| and confirmed analytically by Tagoshi and Sasaki '27]. The 
most recent progress is an extension of the series to O (v 11 ) by Tagoshi, Tanaka and Sasaki |28f. For a test particle in 
circular orbit about a Kerr black hole, the initial progress was again made by Poisson |3(|. The series approximation 
was improved from O (v 3 ) to O ( v 5 ), and subsequently to O (v 8 ) , by Tagoshi, Tanaka, Shibata and Sasaki |3ltl32| . 

Several authors |28ll33l I34L 135. 36] have shown that the convergence of both post-Newtonian approximation and 
black hole perturbation theory is too slow to be useful in constructing accurate templates. More recently, Damour, 
Iyer and Sathyaprakash (hereafter DIS) showed for the case of a test-mass in orbit about a Schwarzschild BH, that 
by using properly defined energy and flux functions that have better analytical properties, combined with Pade 
techniques, it was possible to take the existing series expansion and improve its convergence properties |36| . The new 
approximation in which Pade approximants of new energy and flux functions are used to derive improved templates is 
called P-approximant. Using a fiducially defined "exact" waveform, it was shown that the P-approximant templates 
were both more effectual (i.e. larger overlaps with the exact waveform) and faithful (i.e. larger overlaps with the 
exact waveform and lower biases in the estimation of parameters) than the corresponding post-Newtonian (hereafter 
T-approximant) templates. While in general, more templates are needed for P-approximant templates to cover the 
same volume of parameter space |37j , the extra computational cost is preferred for the increased performance in P- 
approximants. The failure of the PN expansion to converge sufficiently quickly in the case of a test particle orbiting a 
Schwarzschild BH 38] does not bode well for the modelling of even a test particle inspiralling into a Kerr BH. Another 
motivation for this work is that at present there is effort to use effective one body (EOB) waveforms [39j,|40| to detect 
the inspiral and merger signals from two comparable-mass Kerr BHs in GW data. As the EOB waveforms are based 
partially on P-approximants, any development and concretization of the benefits of P-approximant templates will 
boost our confidence in using these improved waveforms as detection templates. 



C. Organisation of the paper 

In this paper we will extend the P-approximant technique to the case of a test particle orbiting a Kerr black 
hole. The reason for focusing on test-mass systems is that we can use the exact numerical fluxes JS3 from black 
hole perturbation theory with which to compare our results and thereby reliably demonstrate the usefulness of the 
technique. We begin in Sec. UTI with a summary of the matched filtering and the signal-to- noise ratio achieved by the 
first generation of GEO, LIGO and VIRGO interferometers for spinning black hole binaries. 

We then go on to discuss in Sec. IIIII the current state-of-the-art in our understanding of the evolution of a test 
particle in orbit around a Kerr black hole. In particular, we shall discuss the time-evolution of the orbital energy and 
gravitational wave flux as a function of the spin of the central black hole at various post-Newtonian orders, and the 
locations of the last stable and unstable circular orbits. We shall see that the post-Newtonian expansion of the flux 
does not show a regular behaviour as we move from low to high orders in the post-Newtonian expansion, becoming 
worse for more rapidly spinning black holes. In order to improve the convergence properties of the flux function, in 
Sec. IIII El we shall introduce a modified form of the flux function and its Pade approximant. We shall demonstrate in 
Sec. II VI the improved behaviour of the Pade approximant, at first graphically and then by showing that the overlaps, 
of the inspiral waveform based on it with the exact waveform, are close to unity. We shall use a number of different 
test systems in our comparison: These will range from systems with dissimilar masses, such as a NS-BH binary, to 
comparable-mass systems, such as a BH-BH binary. To deal with the comparable mass systems, we shall introduce 
in the energy and flux functions a term dependent on the symmetric mass ratio r\ = mi 777,2 /M 2 . While not being 
entirely consistent, because of the fact that no finite mass correction terms are included, it allows us to examine the 
performance of various templates as we move from the test-mass systems to comparable-mass systems. 

The emphasis of the current paper is also on the estimation of parameters. We have carried out a detailed study 
of how good T- and P-approximants are in measuring the parameters of the binary. We shall show in Sec. II VI that 
as a result of not being effectual representations of the exact signal, T-approximants also turn out not to be faithful 
representations either. In other words, the systematic error in the estimation of the parameters caused by the wrong 
phasing of the signal is much larger in the case of T-approximants than in the case of P-approximants. In summary, 
P-approximants are not only effectual but they are faithful too as in the case of non-spinning black hole binaries. 

II. WAVEFORM AND SIGNAL-TO-NOISE RATIO 

In this Section we will discuss the nature of the post-Newtonian waveform from an inpiralling compact binary and 
the response of an antenna to such a signal. We will then use the Fourier transform of the waveform to compute the 
signal-to-noise ratio (SNR) that can be achieved for these signals when they are detected using matched filtering. 

A. The Waveform 

In the transverse-traceless gauge gravitational waves are represented by two polarisation amplitudes h + and h x . 
The response h(t) of an antenna to an incoming signal is expressed as a combination of the two polarisation states 
and the beam pattern functions F + and F x of the antenna as 43] : h(i) = h + F + + F x h x . For a wave from a binary 
of masses mi and 7712 (total mass M = mi + 7712 and mass ratio 77 = mi 777,2 /M 2 ) that is inclined with respect to 
the plane of the sky at an angle i, propagating in the direction (9, </>), (see Ref. (43J for exact definitions), frequency 
F and polarization angle ip with respect to the antenna, the polarisation amplitudes, in the so-called restricted 
post-Newtonian approximation |33j, are 

, / , , v 477M (1 + cos 2 i) ,.. ,.. 

h+(t; M,r,,i) = -L-± '-v%(t) cos cf>{t), (1) 

h x (t; M,r],i) = -^— cos iv 2 F (t) cos (j)(t), (2) 

with vp = (ttMF) 1 ' 3 a velocity parameter, and the beam-pattern functions are 

F+{e,<j>,<ij)) = - (1 + cos 2 6) cos 20 cos 2rp - cos 9 sin 20 sin 2$, (3) 

F x (9,(f),t/j) = - (l + cos 2 9) cos 20 sin 2-tp + cos 9 sin 20 cos 2$. (4) 



Using the above expressions for the beam-pattern functions and gravitational wave amplitudes the response takes the 
form 



h(t) = -^—Cvp(t) cos[(f>(t) + 4> ], 



(5) 



where the amplitude coefficient C and phase 4>o can be assumed to be constant for signals lasting for a short duration 
(say, less than about 30 mins): 



C = 



-\J{\ + cos 2 if Fl + 4(cos if 



F 2 



tan 



2F X cos i 



F + (l 



(6) 

(7) 



Thus, gravitational wave antennas are not able to extract the two polarisations separately and the data analysis 
problem boils down to matching the time- varying phase 4>(t), and to a lesser extent the amplitude vp, of the antenna 
response function. 

Post-Newtonian theory and the quadrupole formula applied to a binary give the relativistic binding energy E(vf) 
per unit mass, and the flux of the waves T(vf) as series expansions in the parameter vp. Once we have the binding 
energy and the flux we can use the energy balance argument, namely that the flux of gravitational waves is completely 
balanced by the negative rate of change of the binding energy [M dE{vp)/dt — —T{vp)], in order to arrive at a 
parametrized equation for the evolution of the phase <p{t) of gravitational waves. Integrating the energy balance 
equation supplemented by lixF — d(j)/dt, one obtains: 



r rof e'(v) 

t(vp) = tret + M / dv—^-, 

rV ™* 3 E'(v) 
dvv —--. 



(j){v F ) = rof + 2 



(8) 
(9) 



where E'(vf) — dE(vp)/dvF, i re f is a reference time and </> rc f is a reference phase of the signal, when the velocity 
parameter is vf — ^rcf- The numerical integration of the above equations is more economical when the following 
differential equations are solved instead, 



dvF 
~~dT 



Hvf) 

ME'{v F y 



d0 
~dt 



2vp 
M ' 



(10) 



The parametric representation, Eqs. © and ©, of the phasing formula </> = </>(£) holds under the assumption 
of 'adiabatic inspiral', i.e., that gravitational radiation damping can be treated as an adiabatic perturbation of a 
circular motion. However, the effective one-body approach [39|, |40j, |4lj, |42| has allowed a treatment of the radiation 
damping to proceed beyond the adiabatic approximation. In order to extract an inspiral signal that may be buried in 
noisy data by the method of matched filtering, we need to employ post-Newtonian accurate representations for the 
two functions E'(vp) and T{vp) that appear in the above phasing formulas. Given an approximant Ea{vf), Fa^f), 
one can define by replacing E(vp) —>■ Ea(vf), F{vf) —■ * J~a(vf) hi Eqs. JSJl and J5J some approximate parametric 
representation, t = ^(vp), <fi = (J)a(vf), and therefore a corresponding approximate template 



h A =h A (t;C, t mi , 0ref, M,r)), 

obtained by replacing vf, in the following w^-parametric representation of the waveform 



h A (v F ) 



4r]M 



Cvp cos 4>a(vf) , 



(11) 



(12) 



by the function of time vf = VA(t) obtained by inverting t = tA(vF)- 

There are several ways of performing this inversion which leads to the different T-approximants, TaylorTl [keep 
the rational polynomial E' [v f) / F {v f) in Eqs. © and Q as is and solve the integrals numerically], TaylorT2 [re- 
expand the rational polynomial as a post-Newtonian series, truncate terms to the appropriate order and solve the 
integrals analytically to get series expansions, t — ^2 k t^vp and <j> — ~^2 k (j>kV F , but solve numerically for cf>(t) from the 
parametric equations] and TaylorTS [do as in TaylorT2 but also invert the series expansions to obtain the phasing as 
an explicit function of time: cf>(t) — ^2 k 0fc(i re f — t)~ k ^ 8 ], as discussed in Ref. |41| . 



It is often convenient to deal with the Fourier representation of the waveform. In the stationar y p hase app roximation 
the Fourier transform, defined as h(f) = f^° h(i) exp(2xift)dt, for positive frequencies reads [43. 144, lia. Eg 

h(f)= I /i(i)exp(27™/£)d£ = ^^-^ e ^ (/) ~ f ^ ( 13 ) 

J-oo d \Jp 

and, since h(t) is real, h(—f) = h*(f). Here tf is the stationary point of the phase ip (i.e., tf is defined by dip/dt\ t} = 0), 
Vf = (7TM/) 1 / 3 , F is the time-derivative of the instantaneous gravitational wave frequency evaluated at the stationary 
point given by, 

F = LlAlLL (14) 

and tp{f) = 2 n ftf ~ 0(^/) is the phase of the Fourier transform in the stationary phase approximation given by 

fl>(f) = 2mft rct - rof + 2 / '° f (v 3 f - v 3 ) ^p- dv. (15) 

Instead of solving the integrals in the above equation it is numerically efficient to solve the following equivalent 
differential equations for the phasing formula in the Fourier domain: 

# o, n dt nM 2 E'(f) 

¥~ 27rf = ' W + 1^TUJ = - (16) 

And in this study we have used the above differential equations in computing the waveforms. For the energy we have 
used the exact analytical formula (see next Section); for the flux we have used the exact numerical flux to define the 
exact waveform and the perturbative expansions, or their re-summed improved versions, to define the approximate 
waveforms. 

B. Signal-to-Noise Ratio 

In order to estimate the signal-to-noise ratio (SNR) we shall assume that the signal is detected using matched 
filtering and that the bank of templates used in matched filtering has a waveform that matches the unknown signal 
very closely. To compute the optimum SNR we need to known the noise spectral density of the detectors. For initial 
LIGO, the one-sided noise power spectral density (PSD) from the design study 0] is given by (46J 

Sh(f) = 9 x 1CT 46 [0.52 + 0.16a;- 4 - 52 + 0.32a; 2 ] Hz" 1 , (17) 

where x = f / fk, and fa = 150 Hz is the "knee-frequency" of the detector. We take the PSD to be infinite below the 
lower frequency cutoff of /i ow = 40 Hz. For the VIRGO detector the PSD is given by 

S h {f) = 3.24 x 10- 46 [I0 3 (25x)~ 5 + 2/x+l + x 2 ] Hz" 1 , (18) 

where fa = 500 Hz and /i ow = 20 Hz. Finally, for the GEO detector the expected noise PSD is J6(| 

S h (f) = 10- 46 x [(3.4a;)- 30 + 34/a; + 20(1 - a; 2 + a; 4 /2)/(l + a; 2 /2)] Hz" 1 (19) 

where fa = 150 Hz and the noise is assumed to be infinite below a lower frequency cutoff of /i ow = 40 Hz. 
Next, let us define the scalar product of two waveforms h and g by 

< a ^ = 2 r ^Tn \MfW(f) + h*(f)g(f)] , (20) 

where the * denotes complex conjugate and h(f), g(f) are the Fourier transforms of h(t), g(t). The (square) of the 
signal-to-noise ratio obtained by matched filtering a signal h with a template g is given by 

2 _fS\ 2 (h\gf 

P= UJ = W (21) 
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FIG. 1: The RMS SNR for equal mass binaries at 100 Mpc for the LIGO (left), GEO (center) and VIRGO (right) antennas. 



If the template perfectly matches the signal then the SNR is simply p 2 

Eq. IJ13JI . and using the Newtonian values for the energy and flux functions, namely E' = —TjVf, T = 32?7 2 u^ /5, and 

F = 96rjv 1 f 1 /(5TrM 2 ), the SNR reduces to: 



(h, h) . For the inspiral signal given in 

.10 

/ 



p(M,r],d;9,(f>,i),i) 



M 5/6 r] 1/2 
dir 2 / 3 



C(8,$,$,i) 



c»t J- -7/3 



df 



1/2 



(22) 



Here F cxli is the frequency where the data analyst chooses to terminate the signal. The cutoff frequency is usually 
chosen to be the frequency of the waves at the last stable orbit of the binary (cf. Sec. lIII A|) which depends on the total 
mass of the two stars and their spin magnitude. Note that the SNRs depend only on the combination M. = Mrj 3 ^ 5 
- the chirp mass - and not on individual masses of the system. For a given total mass the SNR at a fixed distance 
(alternatively the span of the antenna at a fixed SNR) is the highest for binaries of equal masses (when rj = 1/4) 
and reduces by the fraction y/irj for binaries of unequal masses. The SNR decreases inversely as the distance to the 
source and depends not only on the intrinsic parameters M, rj and F cut , but also on the relative orientations of the 
source and the antenna specified by the angles 9, cf>, i, i/i. We can compute the SNR for a source of RMS orientation 
by averaging the amplitude factor C over all the angles. Performing the average over the angles leads to: 



(Fl] 



6,4>,i> 



(F 2 



x/0,<f>,ip 



(C 2 



' i y 0,<f>,tp 



(23) 



For an ideally oriented source, that is for a source that produces the highest SNR, C — 1. Thus, the RMS and ideal 
SNRs, obtained by using C = 2/5 and C = 1, respectively, are given by |6lj 



Prms 



M 5/6 
dn 2 / 3 



2 
15 



at J- 



-7/3 



Sh(f) 



df 



1/2 



Pideal = ^PRMS- 



(24) 



In Fig. [T] we have plotted the RMS SNR for the LIGO, VIRGO and GEO detectors, as a function of the spin 
magnitude and the total mass of the system. Clearly, systems with large spins in which the test mass is in prograde 
orbit produce the highest SNRs since their last stable orbit frequency is higher than those of retrograde orbits and 
they last longer. To capture such signals, however, one must know the phasing of the waves very well since the bodies 
spend a substantial fraction of their time in the detector band in the strongly non-linear regime of the evolution. We 
shall show later in this paper that there could be a significant drop in the SNR, and a corresponding drop in the 
number of candidate events that can be detected, while P-approximants lose little SNR and therefore capture almost 
all events. 



III. GRAVITATIONAL BINDING ENERGY AND FLUX FUNCTIONS 



We can see from Eqs. JHJ and I© that the phase of the gravitational wave depends both on the energy and 
flux functions of the binary system. In the test-mass case, that is when one of the masses is far smaller than the 
other (477 <C 1) so that the dynamics is governed entirely by the Kerr geometry of the massive body, we have an 
exact expression for the energy E{x), but only a series representation for the flux F(x). In the comparable-mass 



case, wherein one cannot neglect the perturbation caused by the companion, both functions are represented only 
by post-Newtonian expansion. The standard approximants for E(v) and !F(v) are simply their successive Taylor 
approximants Ex n and J~T n , respectively. The gravitational wave signal constructed using these Taylor approximants 
are called T-approximant signals or waveforms are simply T-approximants. 

It was shown in Ref [36| that the convergence properties of both the energy and flux functions can be improved by 
using Pade approximation of well-defined energy and flux functions. Damour, Iyer and Sathyaprakash followed a two- 
pronged approach to construct improved energy and flux functions |36| | : Starting from the more basic energy-type and 
flux-type functions, e(v) and l(v), Ref. [3|j constructs Pade-type approximants, say ep n , lp n , of the "basic" functions 
e(v), l(v). One can then compute the required energy and flux functions entering the phasing formula. The successive 
approximants E[ep n ] and T\ep n , lp n ] have better convergence properties than their Taylor counterparts -Ejw [ e T n ] and 

In this study we will restrict ourselves to the test mass approximation and, therefore, employ the exact energy 
function. However, in the same approximation the flux function has been computed analytically only as a Taylor 
series although numerically the flux has been computed exactly. Since our aim is to draw conclusions on how effective 
P-approximants are in the comparable mass case, wherein one has only a Taylor expansion of the flux, we construct 
P-approximants of flux and compare it with the numerical results. The reasons for constructing a new flux function 
are the following: Firstly, it is well known, in the Schwarzschild case, that a simple pole exists at r = 3 M 23] in the 
expression for the flux. We might, therefore, expect that a similar pole exists in the Kerr case. A Taylor approximant 
of the flux function will never produce a pole while the equivalent Pade approximant, which is essentially a rational 
polynomial, will have a pole. Alth oug h it is well known that Pade approximants have the ability to model functions 
which are subject to singularities J43, H3 it is n °t guaranteed that the rational polynomial will reproduce the exact 
pole. Secondly, it is often possible to recover from the divergence |fi2| ] of a Taylor series by using Pade approximation. 
Again, one cannot be certain that the resulting Pade approximant, even when it converges, will be closer to being 
exact than any of the Taylor approximants, but experience with the test mass approximation in the Schwarzschild 
case render some optimism to this expectation. 

We create a Pade approximation of a truncated power series of an analytic function containing n terms by equating 
it to a rational function such that the rational function when expanded as a power series and truncated to order n 
coincides with the original power series. More precisely, Pade approximation can be thought of as an operator Pj^ 
that acts on a polynomial J2k=o a kV k to define a rational function: 

k=0 



*$(£a^U-^ (25) 



Vfc=0 



fc=i 



such that the number N + M + 1 of coefficients in the rational polynomial on the right hand side is the same as the 
number n + 1 of Taylor coefficients on the left hand side. By setting N = M + e with e = 0, 1, we can define two 
types of Pade approximants: These are the super-diagonal, P M , and sub-diagonal, -P^/ +e , approximants. Normally, 
the sub-diagonal approximants are preferred over super-diagonal approximants. This is because when M = N + e the 
right hand side of Eq. (|25|l can be re-expanded as a continued fraction which have the property that as we go to each 
new order of the power series only one new coefficient needs to be calculated. Conversely, with the super-diagonal 
approximants, we would have to re-calculate all the A's and S's in the above equation as we go to higher orders in 
the Taylor expansion. This means that the sub-diagonal Pade approximants are more stable and if we see a trend 
of convergence in the coefficients the addition of a term is not likely to spoil this convergence. We refer the reader 
to Appendix ^ for a more summary of the properties of Pade approximations and how to find the diagonal Pade 
coefficients using continued fraction expansion. 

A. The Orbital Energy. 

In the case of both Schwarzschild and Kerr black holes we have an exact expression for the orbital energy of a test 
particle in a circular orbit around the parent black hole. For a black hole of mass M the energy E in terms of the 
dimensionless magnitude of velocity v = \J 'M/r, r being is the radial coordinate in the Boyer-Lindquist coordinates, 
takes the form 50] 

l-2v 2 + qv 3 , s 

E(v,q)=r,-= \ 26 

VI -Zv 2 + 2qv 3 



where q is a dimensionless spin parameter given in terms of the spin angular momentum J of the black hole by 
q = J/AI 2 = a/M, with a spin angular momentum per unit mass in the Kerr metric. It is actually the derivative of 
the orbital energy that appears in the phasing formula given by: 



E'(v, q) = —VT] 



1 -6v 2 + 8qv 3 -3q 2 v 4 



(l-3v 2 + 2qv 3 ) 



3/2 



(27) 



where a prime denotes a derivative with respect to the velocity parameter v. 

For black holes with spin, the positions of the last stable circular orbit is a function of the spin parameter q. In 
order to find the position of the last stable circular orbit we take E'(v,q) — which gives |50|: 



riso(l)= M 3 + z 2 (q) T \/\? ~ Zl (q)] [3 + Zl (q) + 2 z 2 (q)} 



where 



2\ 3 



zi(q) = l+(l-g 



**(«) = ^3q 2 + z 2 . 



(1 + g) 3 +(l-q) 3 



(28) 



(29) 



(30) 



The + (— ) sign on r^so corresponds to prograde (retrograde) orbits. Now the position of the last unstable circular 
orbit or "photon ring" occurs where the energy function exhibits a pole. Thus, the photon orbit is found by solving 
a cubic equation. Of the three poles, the physically relevant pole is given by J5(J, l51| 



r%{q) = 2M 



1 + cos 



cos l (=F<?) 



(31) 



where once again the different signs define a prograde (retrograde) orbit. In the limit of q — > 1, both r^so and r pr 
move towards the horizon of the black hole, 



H 



m(i± vT 



<r 



(32) 



In the maximally rotating case of q = 1, we cannot distinguish between the last stable orbit and the light ring. This 
means, the greater the spin of the BH, the closer the particle can approach to it before beginning the plunge phase. 
On the other hand, in the limit q — ► — 1, the position of r^so moves outwards from the BH. Thus the particle begins 
its plunge much earlier than in the prograde case. Finally, for a particle in an orbit about a Schwarzschild black hole 
the above equations give the familiar position of r = 6 M for the LSO and r = 3 M for the photon ring. 



B. The Flux Function 



As discussed in the Introduction it has not been possible to derive an exact analytical expression for the flux 
of gravitational waves emitted by a binary system although analytic approximate expressions, and exact numerical 
results, have been computed. In the interesting case of two comparable masses in orbit around each other post- 
Newtonian methods have been used to derive an expression for the flux to quite a high order in the expansion 
parameter. However, since post-Newtonian theory is known to be poorly convergent, especially when the expansion 
parameter approaches unity, it has been suggested to employ in its place an equivalent rational polynomial, or Pade 
approximation, to the (modified) flux function. Since the purpose of this paper is to test the effectiveness of such an 
approximation we need a firm ground on which we can conduct our test. 

In the case of a test mass in an equatorial orbit around a Kerr black hole numerical methods have been used to 
compute the flux to all post-Newtonian orders. This is, of course, valid only in the limit of vanishingly small mass 
of the test body as compared to the central object. However, all the relativistic corrections, including hereditary 
effects, such as the back scattering of gravitational waves off the curved background geometry, would be present in 
this computation. The deformation of the geometry due to the presence of the second body, or the back reaction of 
the waves on the motion of the test body, which would in turn affect the emission process, will not have been included 
in such a calculation. The deformation of the geometry, parametrized by the symmetric mass ratio of the system 
•q, is important only a few orbits before the two objects merge. This should be expected from the fact that both in 
Newtonian and Einsteinian gravity it has been possible to construct an effective one-body formalism to describe the 



dynamics of a binary and the dynamics derived within this formalism is expected to be valid close to the point when 
the two bodies plunge towards each other. Therefore, we expect that the deformation of the system will not bring 
about a major change in the analytic behaviour of the flux and lessons learnt in the test mass approximation will be 
applicable, albeit qualitatively, in the comparable mass case. 

Thus, our strategy is to employ the test mass exact flux computed numerically together with the analytic expression 
for the energy function discussed in the previous Section. These two exact functions, together with the energy balance 
equation, can be used to construct an exact phasing formula for the inspiral signal. We shall also define an approximate 
phasing formula using the corresponding post-Newtonian expansions of the two functions. Our approximations at 
each post-Newtonian order will be one of two types: (1) The standard post-Newtonian or (2) the improved Pade 
approximation. In the next Section we will discuss how these approximations can be further improved and in the last 
Section we will measure the quantitative performance of the two approximations. 

C. Post-Newtonian flux function 

For a test-particle in a circular equatorial orbit the post-Newtonian expansion of the flux function has been calculated 
up to O (v 11 ) in the case of a Schwarzschild BH 28], and to 0(v 8 ) in the case of a Kerr BH J3J, y2(. The general 
form of the flux function in both these cases is given by the expression 



F Tn {x; q) = F N (x) 



Y^ a k (q)x k + \n{x) ^ h(q)x k + O (a 

.fe=0 fc=6 



(33) 



where the expansion is known to order n = 8 and II, in the Kerr and Schwarzschild cases, respectively, q is the 
dimensionless spin parameter introduced earlier and Fn (x) is the dominant Newtonian flux function given by 

F N (x) = B r?x w . (34) 

5 

Here, x is the magnitude of the invariant velocity parameter observed at infinity which is related to the angular 

I/O 

frequency f2 by x = (M $7) . The relation between the parameter x and the local linear speed v in Boyer-Lindquist 
coordinates is given by 

x{v,q)=v[l-qv 3 + q 2 v (i ] 1/3 , (35) 

which reduces, in the Schwarzschild limit, to X = v. Note that this local linear velocity is related to the Boyer-Lindquist 
radial coordinate r by v 2 = M/r. The spin-dependent coefficients, ak(q) 1 and the log-term coefficients, bk, are given 
by 63], 

r , — 1 n — r, — 1247 „ * 11 9 „ 44711 , 33 q 2 n 8191 tt 59<? 

a-Q — i, ai — U, a 2 — — -ggg-, a 3 — 47T — , a 4 — - 9Q72 + -jg-, a 5 — g^2 iff 

_ 6643739519 _ 1712 7 . 16 tt 2 _ 3424 In 2 _ 65 irq , 611 q 2 
" 6 — 69854400 105 ~*~ 3 105 6 + 504 ' 

_ 16285 7T , 162035 q , 65-irq 2 _ 71 q 3 
7 ~ 504 ' 3888 '8 24 ' 

_ 323105549467 , 232597 7 _ 1369 it 2 , 39931 In 2 _ 47385 In 3 _ 359 it q , 22667 q 2 , 17 q 4 
8 3178375200 ~ l ~ 4410 126 " t " 294 1568 14 " + " 4536 ^ 16 ' 

r _ 1712 . _ n h _ 232597 /nc\ 

be 105"' b7 U ' ° 8 ~ 4410 ' ( 6b ) 

where 7 is Euler's number. For graphical purposes, it is more appropriate to deal with the Newton-normalized fluxes 
defined by 

F Tn (x) = F Tn {x)/F N (x). (37) 

D. P-approximant of the flux function 

We will now outline the method of calculation the P-approximant of the GW flux as proposed by DIS j^Q. We 
notice from the form of the series expansion for the flux, Eq. (|33|l . that we begin to encounter logarithmic terms at 
order x 6 and above. In general, series approximations of this form have slow convergence properties. In order to 
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prepare the series representation of the flux for creating the Pade approximation, it is convenient if we factor out the 
logarithmic terms. We can then write Eq. (|33|l as 



FtJx) 



l + ln(^—)JTl k x k 



£ 



C, X 



J Lfc=0 



(38) 



where the new coefficients c^ and Ik are functions of the old coefficients a% and bu- As in Ref. |3(| the log-terms 
have been "normalized" using the value of the velocity parameter at the LSO; this helps in reducing the importance 
of the log-terms. Factoring out the logarithmic terms aids in constructing the rational polynomial, or the Pade 
approximation, of the reminder. Moreover, since we expect the flux to have a pole at the location of the light ring it 
is best to factor out the expected pole so that the reminder has good analytical properties. To this end we create the 
factored flux function, fr n (x) by the operation 



/t„ 



1- 



%pole 



F T 



(39) 



Factoring out the pole also helps to alleviate the problem arising from the absence of the linear term in the PN 
expansion of the flux. (Note that a\ = in both the Schwarzschild and Kerr cases.) We can see from Appendix 1X1 
that in the absence of such a term the continued fraction form of the Pade approximation, the so-called diagonal Pade 
approximant, would lead either to zero or infinite Pade coefficients. The above operation rectifies this problem by 
introducing a linear term into the Taylor series for the flux. If we write the expression in full we obtain 



.Ms) = 



1 + ln 



Uw^T 



\-j j> 



,fe=0 



(40) 



where f = c and f h =c h - c k _Jx po i e , k = l,...,n. 

We can now construct a new flux function by using the Pade approximant of the factored flux given above. Indeed, 
we can construct two variants of the new flux: The first one is what we call the Direct or D-Pade approximant, 
which is obtained by directly starting from the flux function f(v) in Eq. Ij40|l and constructing the equivalent rational 
polynomial. This is the approach followed in DIS. An alternative approach to this is motivated by the fact that in 
the gravitational wave phasing formula the flux appears in the denominator. Thus, instead of constructing the Pade 
approximant of the flux function one could first construct the polynomial expansion of the inverse of the flux function 
and construct the Pade approximant of the resulting polynomial. We call the approximant constructed this way as 
Inverse- or I-Pade approximant because it is obtained from the Taylor expansion of the inverse of the flux function 
f(v) in Eq. I|40|) . Thus, our two improved versions of the flux are defined as follows: The Direct Pade approximant is 
defined by 



fDP„ (v) = 



1 + ln 



xlso 
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fc=6 
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(41) 



where P™ +e is the diagonal or sub-diagonal Pade approximant, n = 2m + e, with e = or 1, depending on whether n 
is even or odd. The Inverse Pade approximant of flux is defined by 



xlso 



f IPn (v)= I -In 
where the coefficients dk in the Taylor expansion are denned by 
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(42) 



(43) 
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One finds d^ by first expanding the RHS in a binomial series and then identifying the coefficients of the various terms 
with those on the LHS. 

Having constructed the rational polynomials equivalent to a given truncated Taylor series of the modified flux 
function we can return to the original flux function that appears in the phasing formula. We shall call the flux 
function so constructed as P- approximant, and not just Pade approximant to remind ourselves that the new flux has 
been obtained by improving the convergence properties in two steps (i.e., definition of a new flux function and the 
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FIG. 2: The convergence of the Newton-normalized T-approximants (left panels) and Inverse P-approximants (right panels) to 
GW flux at post-Newtonian orders 0(x 5 ) to 0(x 8 ), for a test particle in prograde orbit around a Kerr black hole. The solid lines 
represent the exact numerical fluxes with thicker lines corresponding to larger spin magnitudes, q = {0, 0.25, 0.50, 0.75, 0.95}, 
respectively. The correspondingly thickened dashed lines represent the Taylor- and P-approximants. 



construction of its rational polynomial) and not just a direct application of Pade approximation. Thus, wc define the 
Direct P-approximant of flux as 



and Inverse P-approximant as, 



FdpM= 1 



FipAv) 



Xpole 



•Epole 



fop n {v), 



IipM 



(44) 



(45) 



Detailed investigation shows that Inverse P-approximants of flux, namely Fip n , have better convergence properties 
than Direct P-approximants. We shall therefore use only Fjp n in all our investigations in the rest of this paper. 

E. Convergence of the Kerr flux function 

In this Section we will look at the convergence properties of the post-Newtonian series representations and the 
improved P-approximants of the flux function discussed in the last Section. To this end we shall employ all the post- 
Newtonian orders up to which the Taylor expansions have been worked out in the Kerr case. We shall investigate the 
convergence of the approximants in the case of a test particle (i.e., 77 — )■ limit) either co-rotating (i.e., q > 0, prograde 
orbits) or counter-rotating (i.e., q < 0, retrograde orbits), with respect to a central black hole whose spin parameter 
q takes nine values over the full range from (—1, 1). In this Section we shall make a graphical comparison of the 
analytical and numerical fluxes and assess the performance of the flux with respect to variation in the spin-parameter 
values at increasingly higher post-Newtonian orders. We will find that neither the post-Newtonian series nor their 
improved versions fit the numerical flux exactly when the parameters are chosen to be the same for analytical fluxes as 
the numerical flux. However, the P-approximant fluxes agree with the numerical fluxes quite closely when their spin 
values are slightly mis- matched with the true values of the numerical fluxes. The post-Newtonian approximants do 
show this improved behaviour, albeit not to the extent of P-approximants. Thus, we can expect the P-approximants 
to define better signal models than post-Newtonian approximants with regard to both faithfulness and effectualncss. 



ude Orbits, q > 



The numerical fluxes have been computed exactly |52j , but in the test mass approximation, at a set of test values of 
the spin parameter q. For prograde orbits the flux has been computed at q — {0.25, 0.50, 0.75, 0.95}. For the sake of 
completeness, we shall also include in our comparison the numerical flux for a particle in orbit about a Schwarzschild 
BH, i.e. q = 28]. In Fig. [21 we have plotted the Newton-Normalized numerical fluxes alongside analytical fluxes for 
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FIG. 3: The convergence of the Newton-normalized T-approximants (left panels) and Inverse P-approximants (right pan- 
els) to GW flux at post-Newtonian orders 0(x 5 ) to 0(x 8 ), for a test particle in retrograde orbit around a Kerr black 
hole. The solid lines represent the exact numerical fluxes with thicker lines corresponding to larger spin magnitudes, 
q = {0, —0.25, —0.50, —0.75, —0.95}, respectively. The correspondingly thickened dashed lines represent the Taylor- and 
P-approximants. 



both the truncated Taylor approximants (four panels on the left) and (Inverse) P-approximants (four panels on the 
right). For each approximant, we consider four post-Newtonian orders at {x 5 , x e , x 1 ' , x s } and five spin values. The 
solid lines with increasing thickness correspond to the exact numerical fluxes at progressively larger values of the spin 
parameter starting at q = 0, (thinnest line) and ending at q = 0.95 (thickest line). The corresponding dashed lines 
represent the analytical approximant to the GW flux. 

Let us first concentrate on the post-Newtonian fluxes: We notice immediately that at all values of q, and large 
values of x, the 2.5 post-Newtonian approximant 0(x 5 ) is highly divergent. Indeed, this order has been known 
to exhibit improper behaviour in the Schwarzschild case and continues to be so for prograde orbits in Kerr. The 
situation improves at 0(x 6 ), but begins to worsen again at 0(x 7 ), becoming highly divergent at 0(x 8 ). The plots 
serve to illustrate the principal problem with the PN expansion: Going to higher orders of the approximation does 
not necessarily correspond to a better accuracy. 

The sequence of panels on the right-hand-side of Fig. [21 demonstrates that the P-approximant flux has improved 
convergence properties over the post-Newtonian flux for the same range of post-Newtonian orders and spins. The 
extreme divergence which the post-Newtonian flux exhibited at 0(x 5 ) and 0(x 8 ) has been cured. Let us note, 
however, that sometimes the rational polynomial approximation of a Taylor series can produce spurious poles in the 
region of interest 64]. An example of this can be seen at 0(x 7 ) which shows a spurious pole at every value of q. 
The main factor that stands out in the right panels of Fig. [21 is the following: The P-approximant fluxes display the 
same characteristic shape as the numerical fluxes although they don't agree with the exact flux at higher values of 
the parameter x. 

In the case of prograde orbits the P-approximants graphically display a closer convergence to the numerical fluxes: 
While being better than the Taylor approximation, they are still not as good as we would like. However, we must 
remember that what is more important is how the phasing of the waveforms defined by the two approximants match 
with that defined by the exact flux, rather than a simple graphical agreement of the flux. The graphical display is 
useful as it gives us a rough idea on any improvement in performance of a particular model. But as the noise of the 
interferometer weights the inner-product between templates, it may emphasize a particular range of frequency of the 
waveform more than the other. It may turn out that the divergent region of a particular model is outside the detector 
bandwidth for most masses. In such cases the approximant may work well in spite of the divergences or poles. We 
shall return to this question in Sec. IIVI when we consider the matching of waveforms. 



2. Retrograde Orbits, q < 



We now focus on the GW flux from a test-mass counter-rotating with respect to the black hole. We will again use 
a similar set of spin values, i.e. q = {—0.25, —0.50, —0.75, —0.95} , as well as the Schwarzschild flux for comparative 
purposes. In Fig. [21 left panels, we plot the Taylor approximant flux for retrograde orbits. The solid lines in 
the figure correspond to the exact numerical fluxes with thicker lines corresponding to smaller spin values q = 
0, —0.25, —0.50, —0.75, and —0.95, respectively. The correspondingly thickened dashed lines represent the Taylor 
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TABLE I: The "best-fit" value of the spin parameter, q, for T and P-approximant fluxes for q — —0.95 to q = 0.95. For the 
Schwarzschild case we have an absolute error of 0.11 in the case of T-approximant and —0.08 in the case of P-approximant. 



Spin Magnitude 


-0.95 


-0.75 


-0.50 


-0.25 0.25 0.50 0.75 0.95 


q T 

qp 


-0.90 
-0.99 


-0.68 
-0.81 


-0.42 
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series approximation to the GW flux with the same spin magnitude as the numerical flux. The right panels are the 
same except that in place of Taylor approximant fluxes we plot the Inverse P-approximant fluxes. 

The first thing we notice is that just as in the case of prograde motion the (D(x 5 ) and 0(x s ) approximants are again 
highly divergent although at 0(x s ), the convergence of the PN approximation improves as we go to larger negative 
spins. This is not entirely surprising: As we go to larger counter-rotating spins the higher order Taylor expansion 
coefficients become smaller and the non-linear general relativistic corrections are not as important as in the case of 
prograde orbits. Moreover, the position of the LSO moves outwards from the BH as q becomes more negative. At 
q = —0.25, the LSO is at r^so ~ 6.8 M. This corresponds to xlso ~ 0.38. At q = —0.95 the LSO has moved 
outwards to tlso ~ 9M, corresponding to xlso ~ 0.33. If we compare this to the equivalent spins in the prograde 
case we find that for q = 0.25, rzso ~ 5M and xlso ~ 0.44. At q = 0.95, r^so ~ 2M and xlso ~ 0.7. In 
other words, retrograde orbits are weakly bounded orbits in comparison to prograde orbits and experience smaller 
relativistic corrections. Recalling that gravitational wave luminosity ex x 10 , retrograde flux for q = —0.95 at LSO 
is smaller by a factor (0.70/0.38) 10 ~ 500 compared to prograde flux for q — 0.95 at LSO. Consequently, the PN 
approximation seems to work well in the retrograde case. 

In the case of retrograde orbits P-approximants provide, at order 0(x 5 ), only a marginal improvement over T- 
approximants. This is in sharp contrast to the prograde case where the improvement was excellent. At orders 0(x 6 ) 
and 0(x 8 ) are again quite good, with no poles at 0(x 7 ) for most spin values. While there seems to be a less obvious 
benefit in using the P-approximant flux in the retrograde case, it is again important to note that the P-approximants 
show the same general behaviour as the exact flux in the retrograde case too and can therefore be expected to match 
the true general relativistic signal more closely. 

F. Improving the performance of the Taylor- and P-approximants 

In the last two Sections we noted how the P-approximant sequence seems to have a character similar to the sequence 
given by the numerical computations although the two do not agree perfectly. This suggests that if the value of the 
spin parameter used in the T- or P-approximant fluxes is varied relative to the true value of the spin parameter in the 
exact flux, we should be able to obtain a better fit to the numerical fluxes. While we will initially do this graphically, 
we remind the reader that this exercise is purely demonstrative. As the fluxes are Newton-normalized, it is difficult to 
say in what region of the velocity we should attempt to make the best fit. Added to this is the fact that when matching 
templates we should take the weighting of the PSD into account. However, if the approximant fluxes agree graphically 
well with the numerical fluxes then it would mean that we should be able to match the T- or P-approximant flux with 
the numerical result without introducing any new parameters. To examine how best mis-matching spins might work 
we will focus on one of the post-Newtonian orders, say 0(x 6 ), the order at which the T-approximants work well and 
P-approximants have no spurious poles. 

In Fig. 01 left panel, we have plotted the best-fit for the T-approximant flux for various values of the spin magnitude 
from q = —0.95 (top most curve) to q — 0.95 (bottom most curve). For each spin magnitude we have fitted the T- 
approximant with the numerical result as best as possible by varying the spin-magnitude of the approximant. In 
best fitting the numerical flux we incurred errors (defined as a q = 1 — gfittcd/<Zexact) m the spin magnitude of the 
approximant relative to the exact case. These are listed in Table H] Although the errors decrease as we approach the 
extreme prograde flux, we can see from the figure that the fit is not very good for the Taylor approximant for q > 0. 
Schwarzschild case. 

We have plotted the best-fit results for the P-approximant in right panel of Fig. 0] We can observe that while the 
fit is not as good as the PN flux in the extreme retrograde case, it gets superior as we move towards the extreme 
prograde flux. In fact, the P-approximant best-fits work extremely well up to a spin of q — 0.5. It also works well at 
q = 0.75, but only at low values of the parameter x. For q — 0.95 we did not obtain a good fit without introducing a 
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FIG. 4: A best-fit for T-approximant (left) and P-approximant (right) fluxes by varying the spin parameter q at the 0(x 6 ) 
level for q — —0.95 (top most curve) through to q = 0.95 (bottom most curve). Once again, the solid lines represent numerical 
fluxes and the broken lines represent the approximated fluxes. 



new parameter. It is clear that there is an advantage in using the P-approximant flux. However, we will have to wait 
until we calculate the fitting factors between templates to see just how important this advantage is. 

Let us conclude this Section by noting that when approximate fluxes, both T- and P-approximants, best fit the 
exact numerical flux, errors in relative spin magnitude are large at low values of spins, i.e. several 10's of percents for 
\q\ < 0.50, while the errors are low at large spin magnitudes, i.e. a few percents for \q\ > 0.75. We will see that this 
will exactly be the case even when we compute the best overlap of the approximate signal with the exact signal. 

IV. EFFECTUALNESS AND FAITHFULNESS OF T- AND P-APPROXIMANTS 

In Sec. |n] we have seen that the phase of the GW depends on the evolution of the binary's binding energy and 
gravitational wave flux functions. As the energy is known exactly for a test-particle in circular orbit about a Kerr 
black hole it is important to compute the flux to sufficiently high accuracy so that we can match the phasing of the 
exact waves to a good precision. Investigations of the PN approximation showed that while T-approximants have bad 
convergence properties P-approximants improved the convergence of the flux. We also found that we could use the spin 
parameter q as a free parameter in order to obtain a better matching of the flux. By best fitting the spin parameter, we 
were able to graphically achieve a closer fit to the numerical fluxes. Equipped with the new, improved, P-approximant 
we now move onto explore how well a waveform based on it matches the phasing of the exact waves. We shall address 
the performance of the approximants in extracting the exact waveform in two ways: The effectualness of the templates 
measured in terms of the maximum overlap they can achieve with the exact waveform when the parameters of the 
approximant are varied in order to achieve a good match. The faithfulness of the approximant templates measured 
in terms of the systematic errors in the estimation of parameters while detecting exact waveforms. 



A. Overlaps and fitting factor 



We can employ the results of matched filtering, and its geometrical interpretation, to assess how well our approx- 
imant waveforms match the exact waveform. The geometrical theory of signal analysis |54l l5q can be summarized 
as follows: The set of all detector outputs, each lasting for a duration T, and sampled at a rate f s , and consisting of 
N = f s T samples, can be thought of as a linear vector space. Parametrized signals, such as a binary inspiral waveform 
that depends on the two masses of the component stars and their spins, are also vectors in the vector space of all 
detector outputs. However, the set of all waveforms do not form a vector space although they do form a manifold 
with the parameters serving as coordinators. 

Matched filtering technique, which is the most effective technique to capture signals of known phase evolution, can 
be used to define a metric on such a manifold. The metric and the associated scalar product allows us to compute 
the distance between any two signal vectors that are parametrized in the same way. The scalar product is the same 
as that introduced in Eq. I|20|) . If we have two different families of waveforms, for example T- and P-approximants, 
parametrized similarly then each family lives on a distinct manifold. The distance between two vectors with exactly 
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FIG. 5: The maximized prograde overlaps for T-approximant (left) and P-approximant (right) templates at the x s approxima- 
tion. Each system consists of a 1.4 Mq NS inspiralling into a central BH of mass ranging from 10-50 M@. The figure covers 
prograde Kerr from q = 0.5 onwards, as fitting factors of > 0.99 are achieved for all lower spins. 



the same parameters but belonging to different families is, in general, not zero. Although the scalar product can be 
computed between any two vectors in the space it is only in the case of two similarly parameterized signal vectors does 
the scalar product represent the distance between the vectors. This is because the coordinates, that is the parameters, 
don't have any meaning for vectors that don't live on the manifold. By calculating the scalar product between two 
normalized vectors, that is vectors whose scalar product with themselves is unity, we can see how well they match. 
For two normalized waveforms, or signal vectors, the scalar product returns the cosine of the angle between them and 
is is normally referred to as the overlap, denoted by O. Given two waveforms h and g, not necessarily belonging to 
the same family of approximants, their overlap is defined as 



O 



(h\g) 



V(h\h)(g\g) 



(46) 



where, the inner-product between two real functions h(t) and g(t) is defined by Eq. H2U|) . For detection of signals what 
is more important is the fitting factor FF : As each template is a function of the intrinsic parameters A M , the fitting 
factor is defined as the maximum overlap obtained by varying the parameters of the template (or the approximate 
waveform) relative to the exact waveform: 



FF = max O (Xn . 



(47) 



If two waveforms are a perfect match then their overlap is unity. As the waveforms begin to differ, their overlap 
differs from 1 and the value of the overlap is a measure of how similar the two waveforms are and how useful one 
waveform is in extracting the other waveform buried in noise. Now the GW from a binary is determined by a number 
of parameters: The two masses, spins, eccentricity and the orientation of the angular momentum at some initial time. 
These parameters essentially determine the dynamics of the system and are called intrinsic parameters 55]. The 
observation of such a system introduces five other parameters: These are the instant io of coalescence (or the time 
at which the orbital frequency reaches a certain value) and the phase 0o of the system at that instant, the angular 
position of the system in the sky, the distance of the binary (or, equivalcntly, the amplitude of gravitational waves at 
Earth) and the polarisation angle (or, equivalently, the ratio of /i+ and h x )- These are called the extrinsic parameters 
|55J| and do not play any role in the dynamics of the system. For a signal that lasts only for a few mins, as would be the 
case for systems expected to be observed in a ground-based interferometer, the direction and the polarisation angle 
cannot be measured in a single interferometer and the distance is only an amplitude parameter. Thus the templates 
are only needed for parameters to an d </>o- I n the case of circular equatorial orbits of a test mass around a central 
black hole, the wave is parameterized by the (io, 3>o> M, n, q). 

In order to improve the detection probability, we would like to maximize over as many of the parameters as possible. 
Maximization over t is achieved by simply computing the correlation of the template with the data in the frequency 
domain followed by the inverse Fourier transform. This yields the correlation of the signal with the data for all 
time-lags. However, one has to worry about circulation correlation effects, and the consequent corruption of the 
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FIG. 6: The 4-PN T-approximant (left) and P-approximant (right) analytical flux function for spins of q = 0.5,0.75 and 0.95 
against the numerical fluxes for the same spin values. We can see that, beyond q = 0.5, the T-approximant flux goes to zero 
before the LSO is reached. This happens sooner, the higher the spin value we go to. This is not a problem for the P-approximant 
flux. 



correlation when the template gets split at the boundaries of the sample sets, but these are easily handled by padding 
the template with sufficiently large number of zeroes before performing the Fourier transform. 

It was pointed out by Schutz |l2l l57| that the maximum of the correlation C of data with a template over 4>o can 
be computed using just two templates - an in-phase and a quadrature-phase template. In other words, generate two 
orthonormalizcd templates with phase $o — and <I>o = tt/2. The maximized overlap is then given by 



maxC = 

*0 



Gl 



C ir/2' 



(48) 



where C = (h A (3> ) | h x ) , C = (h A (*„ = 0) | h x ) and C^ /2 = (h A ($ - tt/2) | h x ). Here, h x denotes some 
normalized "exact" waveform. However, numerically, this procedure does not fully orthogonalize the two templates. 
For this reason we specifically orthonormalize the two waveforms. Using two un-normalized waveforms h — h (<5> = 0) 
and hn/2 = h ($0 = 7r /2), we generate two orthonormalized waveforms according to 



H = ^- 

\ho\ 



H^/2 — — 



K/2 - (K/2 \H ) H 



h w /2 - \hn/2 Wo) H \ 



(49) 



Alternatively, one could define the quadrature-phase template using H^/2 = iHg, which is explicitly orthogonal to 
the in-phase template. The (square of the) maximum of the overlap over $0 is given by the sum-of-squares of the 
correlation with the in-phase and quadrature-phase templates: 



max (C) = J{Hv\h x ) 2 + (H„ /2 \h x Y 

3*0 



(50) 



Thus, the correlation can be easily maximized over the unknown time-of-arrival and the constant phase of the tem- 
plates. However, we can also maximize over all other parameters using a grid of templates. 

In this paper, as we are working in the test-mass approximation, we assume that our system is composed of objects 
with a small mass ratio. For concreteness we assume that the system comprises a 1.4M© NS inspiralling into a central 
BH of varying spin magnitude and mass. Beginning with a central BH of IOMq, we work our way upward to a 50AfQ 
BH. In terms of the symmetric mass ratio, this gives us a range of test systems from r\ = 0.027 for the heaviest binary 
in our list [i.e. (50-1. 4)M Q system] to r\ — 0.11 [i.e., (10-1. 4)M Q system]. We will also look at the limiting case of a 
10-10M© equal-mass system. Strictly speaking, the formulas for energy and flux functions used in this study are not 
applicable to the comparable mass case since we have neglected the finite mass correction terms in these quantities. 
However, the results of such a study should give us an indication of how strong are the relativistic corrections, as 
opposed to finite-mass corrections, in the case of comparable masses. In all cases we vary the spin magnitude of the 
central BH from q = -0.95 to 0.95. 

In the next two sub-sections, we will focus respectively on prograde and retrograde systems. As well as presenting 
the results for the fitting factors achieved, we will also focus on the problem of parameter extraction using templates 
constructed from the Taylor- and P-approximants. Since the main focus in this study is test-mass approximation we 
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FIG. 7: The percentage error in the estimation of the chirp-mass, M c , for T-approximant (left) and P-approximant (right) 
templates at the x s approximation. The values of M c correspond to a 1.4 Mq NS inspiralling into a central BH of mass ranging 
from 10-50 Mq . The figure covers retrograde Kerr, Schwarzschild and prograde Kerr systems. 



shall be interested in only one of the mass parameters and consider errors in estimating the chirp-mass Ai — mrf' 5 in 
addition to the spin magnitude of the central BH. In view of economy we shall only present the results for the highest 
PN order available, namely 0(x & ) order. In all cases our fiducial exact signal h x will be that obtained by using the 
exact expression for the energy in Eq. (|26J) and the exact numerical fluxes generated using black hole perturbation 
theory [52| , and the template will be the approximate waveform constructed using the exact expression for the energy, 
as before, and an approximate expression for the flux, either the T-approximant flux, the corresponding template 
denoted /it, or the P-approximant flux, the corresponding template denoted hp. 



B. Prograde Orbits — Effectualness 



In the first instance, we maximized the overlap only over to and <I>o, while holding all other parameters of the signal 
and the templates the same. The overlaps so obtained were not good enough. In the case of T-approximant templates 
the overlaps at q = 0.25 were in the range 0.6 < FF < 0.8, depending on the PN order. By q = 0.95 they had fallen 
to 0.3 < FF < 0.5. The P-approximants actually fared much worse in this case. At q = 0.25 the fitting factors lay 
between 0.25 < FF < 0.65, while at q — 0.95, none of the fitting factors achieved was > 0.4. 

It is clear that maximizing over extrinsic parameters is not adequate: When we maximize over all other parameters, 
overlaps improve significantly as compared to the unmaximized overlaps. For T-approximants - Fig. \5\ left panel - 
maximizing over all parameters gives fitting factors of FF > 0.98 at all mass ranges for the test-mass systems up 
to a spin of q = 0.75. Between 0.75 and q = 0.95, the T-approximant templates begin to perform badly and the 
fitting factors drop to FF ~ 0.82. For the equal-mass case - Fig. [5]- the templates once again achieve fitting factors 
of FF > 0.99 up to a spin of q = 0.75, but fall off at higher spin magnitudes achieving a fitting factor of ~ 0.98 
at q = 0.95. We should point out that these results do not properly convey just how bad the 4-PN T-approximant 
template actually performs. In the left hand panel of FigEl we plot the PN approximation for the flux function 
against the numerical fluxes at spins of q — 0.5, 0.75 and 0.95. We can see that the flux function at q = 0.75 and 0.95 
become negative long before the LSO is reached. This means that as we go to higher and higher spins, we can model 
less and less of the waveform. For q < 0.6 we can model the waveform up to the LSO. However, for q > 0.6 this we 
have to stop the waveform generation at a cutoff velocity of x cu t ~ 0.5. At q — 0.75 we model approximately 80% of 
the waveform, while at q = 0.95, this reduces to 44%. Thus, T-approximants have to be prematurely terminated well 
before the LSO is reached. In the process of maximization of the overlap we systematically find that exact waveforms 
corresponding to a given total mass are maximized by T-approximant templates of significantly smaller masses. The 
reason for this is that smaller mass templates will have a longer duration (and therefore make up for the loss in overlap 
due to premature termination of the waveform) and thereby achieve a better overlap with the exact waveform. 

Therefore, the fitting factors beyond q = 0.75 completely overestimate the performance of the T-approximant 
template. 

If we now focus on the right hand panel of FiglHl we see the true power of the P-approximant templates. The 4-PN 
template suffers none of the divergences that effect the T-approximants. We therefore generate all templates up to 
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FIG. 8: The percentage error in the estimation of the spin parameter, q, for T-approximant (left) and P-approximant (right) 
templates at the x 8 approximation. The values of M c correspond to a 1.4 Mq NS inspiralling into a central BH of mass ranging 
from 10-50 Mq. The figure covers retrograde and prograde Kerr systems. 



the LSO or 2 kHz, whichever is reached first. We can see from the right panel of Fig. [SJthat the P-approximant 
templates achieve fitting factors of > 0.99 at all spin values. For the equal-mass case - Fig- El — the P-approximant 
templates achieve fitting factors of FF > 0.995 at all mass and spin levels. This demonstrates that in the case of 
prograde orbits, the P-approximant templates are clearly more robust, even at high spin magnitudes of q — 0.95. 



C. Prograde Orbits — Faithfulness 



Gravitational wave observations of the ins pire d signal are expected to result in very precise measurements of the 
signal parameters. Indeed, it has been shown 53] that the chirp mass can be measured to a relative accuracy of 0.02% 
and 0.16% respectively for systems comprising (10, 1.4)M© and (10, 10)M Q objects when the correlations between 
the spin magnitude and masses is neglected and to an accuracy of 0.19% and 1.42% for the same systems when the 
correlations are included. Such accurate estimation of parameters will aid in high precision tests of general relativity. 
When we use approximate templates the question naturally arises if the bias in the estimation of parameters 65] is 
smaller than the accuracy with which the parameters can be determined. This is what we shall consider next. 

While the T-approximant templates performed well up to q = 0.75, a clear indication of how good they truly are 
comes from the bias in parameter estimation. This tests the faithfulness of the templates. It should be pointed out 
at the outset that there is a large covariance between the chirp mass and the spin magnitude and a small error in 
the estimation of the chirp mass is compensated by a large error in the estimation of the spin magnitude. In Fig. [7| 
we have plotted the percentage bias in the estimation of the chirp-mass for both T- (left panel) and P-approximant 
(right panel) templates. From Fig.0 left panel, it is clear that for T-approximants the bias in M. varies between - 
5%. Comparing the left and right panels of Fig. we find that the P-approximant templates (right panels) are more 
faithful with the percentage bias being less than 1% in general. 

In Fig.[S]we present the percentage bias in the estimation of the spin magnitude of the central BH. We find that for 
both approximants we have to endure a large bias in q. Once again the P-approximant templates are more faithful. 
For T-approximants the bias varies between 5 and 90%, while for P-approximants it lies between 2 and 10%. From 
the middle and bottom Sections of Fig. EH it is clear that even in the equal mass case the bias in the estimation of 
both parameters is greater for T-approximant templates than P-approximants. 

It is therefore clear that in the case of prograde orbits, P-approximant templates must be used in any detection 
strategy. Not only do they achieve higher fitting factors at all spin and mass ranges compared to the T-approximant 
templates, but they are also more faithful. Let us finally note that the bias in the estimation of spin is larger at low 
spin values, which is exactly what we found by best fitting the approximate fluxes with the exact fluxes by using q as 
a free parameter. 
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D. Retrograde Orbits — Effectualness 

In the case of retrograde motion we carried out the same range of tests as in the prograde case. We found that in 
the case of maximization over to and $o the overlaps are not as bad. It seems that all templates perform reasonably 
well for all spin values. This is probably due to the fact that the PN coefficients at higher orders are smaller when 
q < than when q > 0. Additionally, retrograde orbits are less bound and, therefore, less relativistic, and again higher 
order PN corrections do not play as significant a role for retrograde orbits as they do for prograde orbits. At a spin 
value of q — —0.25, the T-approximant templates achieved fitting factors in the range 0.67 < FF < 0.81. This falls 
off to a mean value of ~ 0.54 at q = —0.95. The P-approximant templates also perform better in the retrograde case 
achieving fitting factors in the range 0.47 < FF < 0.80 at q = -0.25, and between 0.54 < FF < 0.84 at q = -0.95. 

In addition to extrinsic parameters when we maximize over all the intrinsic parameters the fitting factors improve. 
In fact, in the case of retrograde orbits both templates achieve fitting factors of FF > 0.99 regardless of the spin 
magnitude and the chirp mass. There is no surprise here: The retrograde waveforms are still well within the adiabatic 
regime and are, therefore, modelled well by both templates. From a purely detection point of view, unlike the prograde 
case, there is no obvious benefit from employing P-approximant templates. For the equal mass system too there is 
not much difference in the fitting factors of the two families of templates with the exact signal. However, the fitting 
factor for T-approximants is well below than the corresponding P-approximant template for spin magnitudes q > 0.5. 

E. Retrograde Orbits — Faithfulness 

The benefit of using P-approximant templates for retrograde motion is only observed when we consider parameter 
extraction. Referring to the bottom portions of the two panels in Fig. [7J we note that the T-approximants perform 
well at all spins for 3.0 > M. < 4.5, with a bias of less than 1% in the estimation of M.. Beyond this, there is in general 
a bias of > 2%. As we approach the extreme retrograde case the bias rises to as much as 55%. The P-approximants 
perform in a similar manner. The bias in the region 3 > M. < 4.5 is again in general less than 1%. The error does 
again increase as we head towards the extreme test-mass range, but in this case it reaches a maximum value of 8% as 
opposed to the 55% seen in the case of the T-approximants. 

The pattern observed in the retrograde case is similar to the prograde case with the errors in the T-approximants 
greater than those in the P-approximants. The error in estimating q reaches a maximum of 80% for T-approximants 
and 25% for P-approximants. There is a slight difference, however, between the two cases: We observe that there is an 
asymmetry in the distribution of contour lines. The T-approximant templates suffer a higher error over a larger s pan 
of spins in the retrograde case than in the prograde case. This is consistent with the results presented in Ref. [3lll32]| 
where it was shown that the PN approximation waveforms perform worse in the retrograde case. We should, therefore, 
not be surprised that the errors are greater for retrograde motion. Now, referring once more to the equal mass case, 
the P-approximant templates outperform their T-approximant counterparts. Even with this, we must again conclude 
that while on the surface there is no clear case for using P-approximant templates in searching for retrograde motion 
systems, they should be used because of the lower bias incurred in the estimation of parameters. 

F. Schwarzschild Orbits. 

While we have presented the maximized fitting factors and percentage errors in the estimation of the chirp- mass for 
the Schwarzschild case, we have not made any remarks regarding the maximized spin value. In Table ITT1 we present 
a comparison of the maximized spin for both T and P-approximant templates for a selection of systems. We can see 
from the table that the best overlap is never achieved using the Schwarzschild value for the spin. For T-approximants 
the best fitting factors are all obtained by using retrograde values for the spin. For P-approximants, it varies between 
prograde and retrograde values. In general, regardless of whether the maximized spin is a prograde or retrograde 
value, we can see that the absolute error in the spin value is always smaller for P-approximant templates. 

V. CONCLUSIONS 

We have applied P-approximant templates to the case of a NS orbiting Kerr BHs of varying mass and spin. When 
graphically matching the numerical fluxes we found that the P-approximant fluxes gave a better fit as compared to the 
PN fluxes, especially when we varied the value of the spin parameter and best fitted the exact flux. However, the true 
test of how well a template performs is reflected by the fitting factors achieved. Using a signal waveform constructed 
from the exact expression for the orbital energy and numerical fluxes from black hole perturbation theory, we were 
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FIG. 9: The fitting factors for a 10-10 Mq binary - without the finite mass correction terms - for T and P-approximant 
templates (top) . The percentage error in the estimation of the chirp- mass for T and P-approximant templates (middle) . The 
percentage error in the estimation of the spin of the central black hole for T and P-approximant templates (bottom). 



TABLE II: Fitting factors for a variety of Schwarzschild systems of various masses. The maximized values of mi, m-z and q are 



presented in brackets underneath. 
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able to compare the performance of T and P-approximant templates. In the case of retrograde, Schwarzschild and 
prograde orbits, not only did the P-approximants gave better and more reliable fitting factors, they also gave smaller 
biases in the estimation of parameters. We also saw the true power of the P-approximant templates in that we were 
able to generate templates right up to the LSO . This is something that was not possible with the T-approximant 
templates due to the approximation for the flux function becoming negative before the LSO is reached. This restricted 
just how much of the numerical waveform we could model. While not being completely correct due to the fact that we 
omitted the finite-mass correction terms, we also saw that the P-approximant templates gave the best performance 
when in the equal-mass case. 

The bias in the estimation of parameters is particularly poor in the case of T-approximants especially for high-mass 
systems. Since high-mass systems can be seen to a greater distance, and therefore likely to be the first sources to be 
detected, it is important to keep in mind that the bias in the estimation of parameters could be as large as the random 
errors associated with the measurement. Such biases could seriously affect the test of general relativity J53| that are 
envisaged to be carried out with the aid of high-mass systems. It is clear that for the type of systems examined in this 
paper, namely equatorial test-mass circular orbits in Kerr, P-approximant templates are to be preferred over their 
PN counterparts in both detection and measurement. It remains to be seen if P-approximants (or their improved 
versions) prove to be good enough when considering generic orbits in Kerr and for comparable mass systems with 
spinning black holes on generic orbits. 
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APPENDIX A: PADE APPROXIMATION 

Let f{x) be an analytic function whose Taylor expansion in x, denoted T n (x), is T n {x) = X]fc=o a nX n ■ 
We construct the Pade Approximation Pj^(x) of T n (x) by equating the polynomial expansion with a rational 
function to the same order according to 

T n {x)=P£{x) = %I@r, (Al) 



where 
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such that n + 1 = N + M + 1 . The benefit of using Pade Approximants over Taylor Approximants is that the error in 
the approximation after L terms reduces by a factor of 2 L 49]. If we let M = N + e, where e = or 1, we construct 
the "diagonal" and "sub-diagonal" Pade Approximants respectively. The advantage of this is that we can cast the 
rational function in the form of a continued fraction, i.e. 

-P/V+e = CTW ' (A3) 
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If we use the rational function interpretation, Eq. IA2I each time we go to a higher order n of approximation, we have 
to recalculate each of the coefficients, b n and d n . However, by using the continued fraction formalism, we only have 
to calculate the newest coefficient, c n . As an example, consider the Taylor expansion 

T n (v) = a + ai (v). (A4) 

By matching this to a Pade Approximation of the same order we get 

a + aiv = — . (A5) 

1 +ClV 



Solving for the Pade coefficients, c„, gives us 



c = a , ci = . (A6) 

a 



If we now go to the next order, i.e. 



a + a 1 v + a 2 v 2 = ^ — , (A7) 



1 + c 2 v 



and solve for the coefficients c n once more, we obtain 

c = a Q , c\ = -— , c 2 = ■ -— + — . (A8) 
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The example shows that as we go to higher and higher orders in v, we only have to calculate one new coefficient as 
all the lower order coefficients remain constant. Consequently, there is a sense of stability in the behaviour of the 
continued fraction coefficients and convergence is easier to test. 

There are, unfortunately, problems with Pade approximation that one should bear in mind: Firstly, there are regions 
where Pade approximants diverge while Taylor Approximants converge. The second, and perhaps most important 
feature is that when Pade approximation does indeed converge, we cannot be sure that it converges to the exact 
function that we had in mind |48| . It should therefore always be kept in mind that there are inherent difficulties 
in using Pade Approximants. Finally, Pade approximants sometimes exhibit spurious poles in the region of interest. 
However, as one will be able to test a priori the existence of poles one can avoid the use of such approximants. 



